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ABSTRACT 

In this paper, gravothermal oscillations are investigated in two-component clusters 
with a range of different stellar mass ratios and total component mass ratios. The 
critical number of stars at which gravothermal oscillations first appeared is found 
using a gas code. The nature of the oscillations is investigated and it is shown that the 
oscillations can be understood by focusing on the behaviour of the heavier component, 
because of mass segregation. It is argued that, during each oscillation, the re-collapse 
of the cluster begins at larger radii while the core is still expanding. This re-collapse 
can halt and reverse a gravothermally driven expansion. This material outside the 
core contracts because it is losing energy both to the cool expanding core and to the 
material at larger radii. The core collapse times for each model are also found and 
discussed. For an appropriately chosen case, direct iV-body runs were carried out, 
in order to check the results obtained from the gas model, including evidence of the 
gravothermal nature of the oscillations and the temperature inversion that drives the 
expansion. 

Key words: globular clusters: general; methods: numerical; methods: iV-body sim- 
ulations. 



1 INTRODUCTION 

Gravothermal oscillations are one of the most interesting 
phenomena which may arise in the post-collapse evolution 
of a star cluster. The inner regions of a post collapse cluster 
are approximately isothermal and are subject to a similar in- 
stability as the one found in an isothermal sphere in a spher- 



ical container, as studied by Antonov ( 1962 1 and Lynden- 



Bell & Wood (19681. Gravothermal oscillations, which are 



thought to be a manifestation of this instability, were discov- 



ered by Bettwieser & Sugimoto (1984) whilst studying the 



post-collapse evolution of star clusters using a gas model. 
For a gas model of a one-component cluster it was found 
that gravothermal oscillations first appear when the num- 
ber of stars N is greater than 7000 ( |Goodman||1987[ ). This 
value of N has also been found with Fokker-Planck calcula- 



tions (Cohn et al 19891 and by direct A-body simulations 
(Makino 1996). However, in a multi-component cluster the 
situation is more complicated. The presence of different mass 
components introduces different dynamical processes to the 
system such as mass stratification. Multi-component sys- 
tems try to achieve kinetic energy equipartition between the 
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components, which causes the heavier stars to move more 
slowly and sink towards the centre. This can lead to the 
Spitzer instability (Spitzer 1987) in which the heavier stars 
continuously lose energy to the lighter stars without ever be- 



ing able to reach equipartition. Murphy et al ( 1990 1 found 



that the post-collapse evolution for multi-component mod- 
els was stable to much higher values of N than in the case 
of the one-component system and that the value of A^ at 
which gravothermal oscillations appeared varied with differ- 
ent mass functions. 

In order to gain a deeper understanding of gravothermal 
oscillations, it is desirable to work with simpler models in 
which some of the effects which are present in real star clus- 
ters are ignored or simplified. For example, real star clusters 
have a range of stellar masses present, but in the current pa- 
per, the stellar masses are limited to two. Gaseous models 



are often used in this kind of research (Bettwieser & Sugi- 



|moto|1984(|Goodman|1987||Heggie fc Aarseth|1992|) because 
they are computationally efficient. |Kim, Lee fc Goodman| 
| |1998[ ) have already completed research in this area using 
Fokker-Planck models. However, their research was limited 
to mostly Spitzer stable models and only a small range of 
stellar mass ratios. The study in the present paper looks 
at the more general Spitzer unstable models using various 
stellar mass and total mass ratios. 
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There is also evidence of gravothermal oscillations in 



real star clusters. Giersz & Heggic (2009) modelled the clus- 



ter NGC 6397 using Monte Carlo models and found fluctu- 
ation in the core radius. Their timescale suggests that they 
are gravothermal. Subsequently, they confirmed these fluctu- 
ations using direct iV-body methods with initial conditions 



generated from the Monte Carlo model (Heggie & Giersz 



2009) 



Two-component clusters may seem very unrealistic but 
there is reason to believe that they may be a good approxi- 



mation to multi-component systems. Kim & Lee ( 1997 1 were 



able to find good approximate matches for half-mass radius 
Th, central velocity dispersion v c , core density p c and core 
collapse time t cc between two-component models and eleven- 
component models which were designed to approximate a 
power law IMF. Also see Kim, Lee & Goodman (19981 for 



a discussion of the realism of two-component models. 

This paper is structured as follows. In Section 2, we de- 
scribe the models which are used. This is followed by Section 
3, in which the results concerning gravothermal oscillations 
are given. Section 4 is concerned with the results of the core 
collapse times. In Section 5, the results of iV-body simula- 
tions are given. Finally Section 6 consists of the conclusions 
and a discussion. 



2 MODELS 
2.1 Gas model 

2.1.1 Basic equations and Notation 

In our model, we ignore primordial binaries and stellar evo- 
lution, and assume that the energy generating mechanism 
is the formation of binary stars in three body encounters 
and subsequent encounters of binaries with single stars. In 
a one-component model the rate of energy generation per 
unit mass is approximately 

or G 5 m 5 n 2 
e = 85 =— (1) 



Table 1. Notation (the subscript i corresponds to the i th com- 
ponent, i = 2 refers to the more massive component ) 



radius 
mass density 
one dimensional velocity dispersion 
stellar mass 
total mass (within radius r) 
Constant (see text) 
energy flux 
number of stars 
coulomb logarithm (A = 0.02./V) 

Lagrangian derivative (at fixed M) 



in 
M 
C 
L 
N 
In A 
D 

J 
dr 



radial derivative (at fixed t) 



where i = 1,2. This model in turn is ultimately inspired 
by the one-component model of Lynden-Bcll & Eg gleton| 
(19801). 



The meaning of the symbols can be found in Table [T] 
The major difference between the above equations and those 
for the one-component model is the last term of equation [5] 
which involves the exchange of kinetic energy between the 
two components. See Spitzer (1987, p. 39) for information on 
this term. As the heavier component dominates in the core 
of the cluster, it is assumed that all of the energy is that 
generated from the second component. Hence the Kronecker 
delta #2,i in the last equation. There are two constants in 
the gas code which can be adjusted: C and the coefficient 
A of N in A = XN. The value of A = 0.02 was used as it 
was found to provide a good fit for multi-component models 
(Giersz fc Heggie ||1996[). The value of C used was 0.104 



(Hcg gie~fe Ramamani ||1989 |. This value of C results from 



the comparison of core collapse between gas and Fokker- 
Planck models of single component systems and it is not 
clear if it applies accurately to post-collapse two-component 
models. 



| |Heggie fc Hut|2003l l, where m is the stellar mass, n is the 
number density, a c is the one dimensional velocity dispersion 
of the core and G is the gravitational constant. [Goodman] 
( |1987[ ), whose results on the 1-component model we shall 
occasionally refer to, used a similar formula, with a coeffi- 
cient which is, in effect, in the range 140-170 (depending on 
the value of N). 



The equations of the two-component gas model ( Heggic 
fc Aarseth|1992 | are given below: 

dM, 



c)r 



dr 

doi 
dr 

dU 
dr 



— 4-Kpir 



G(Mi + M 2 ) 



12nCmiPir 2 ]nA 



(2) 
(3) 
(4) 



-47rr 2 p. 



+4(2tt)5G 2 lnA^- 



P3- 



Pi 

(m 3 - 



+ Si. 2 e 



(5) 



_2 
i&3- 



2.1.2 The role of N in the gas code 

This paper places emphasis on the role of N in evolution, 
but it is not clear what role N plays in equations Q - ((HI). 
For fixed structure (i.e. Pi(r), etc), N appears explicitly in 
A (where its role is rather insignificant), and in the individ- 
ual masses mt. These appear in equations Q and |5j. In a 
system with fixed structure, equation Q shows 



Li oc m In A oc 



InAJV 



reflecting the fact that the flux L is caused by two-body 



relaxation, and its time scale is proportional to 



In 



equation |5j) N plays a similar role in the last term on the 
right, which governs the approach to equipartition. It also 
appears implicitly through e, because of the m dependence 
in equation |TJ. For a system of given structure, its contri- 
bution to L in equation |5| is proportional to A^ -3 (as we 
are assuming that p = mn is fixed and so e in equation |T]) 
is proportional to m 3 ). It would seem as though this term 
is insignificant for large N. In practice, however, the sys- 
tem compensates by increasing the central density so that e 
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Table 2. Critical value of N (N crit ) in units of 10 4 



Time (ii,Th) 



Unstable postcollapse evolution 
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Figure 1. Logarithm of the central density vs time (in units of 
the initial value of t r h) for a two-component gas model, = 2, 

f^- = 1, top: N = 1.5x10* (stable), bottom: 2.5x 10 4 (unstable). 
For initial conditions see Section l3.ll 



plays a comparable role to the relaxation terms (see Section 

poll . 



2.2 Direct A-body 

Direct iV-body simulations were conducted using the 
NBODY6 code (Aarseth, 2003) enabled for use with Graph- 
ical Processing Units (GPU). NBODY6 has a range of fea- 
tures and options such as individual time steps which make 
it an excellent direct iV-body code. NBODY6 is written 
in FORTRAN and is publicly available for download from 
www.ast.cam.ac.uk/~ sverre/web/pages/nbody.htm . 



3 CRITICAL VALUE OF N 

If the value of N is not too large, then, after core collapse, 
the cluster expands at a steady rate (Fig. [I] top). However, 
at a larger value of N the central density (p c ) was found 
to oscillate (Fig. [T] bottom). |Goodman| ( |1987| ) showed that 
for one-component models the steady expansion is unstable 
for large values of N and found that the value at which 
oscillations first appeared is N — 7000. In this present paper, 
the case of two-component models is investigated. 



Figure 2. Contours of log 1() (iV cr ;t) 



3.1 Results of the gas code 

In all cases, the initial conditions used were Plummer mod- 
els (|Phimmer||1911| |Heggie fc Hut| |2003[ ). The initial ve- 



locity dispersions of both components were equal and the 
initial ratio of density of each component was equal at all 
locations. The initial conditions were constructed with dif- 
ferent stellar mass ratios — = 2, 3, 4, 5, 10, 20, 50 and for 
each of these mass ratios, a model with total mass ratios 
Iff = 0.1,0.2,0.3,0.4,0.5 and 1 was constructed. A python 
script was used to run the gas model code over a range of 
values of N for each of the pairs of mass ratios. Each run 
terminated when the time value reached 30 initial relaxation 
times (t i}r h)-The value of the central density was checked for 
an increase in value of 5 percent or more in any interval over 
the time period between 20t i>r h and 30ii ir h- If an increase 
was found, the run was deemed to be unstable and the range 
of N was refined. This process continued until the critical 
value of N (N crit ) at which oscillations first appeared was 
determined (correct to ten percent). The values of N cr u were 
also visually confirmed from the output of the gas code. The 
obtained values of N crit in units of 10 4 are given in Table [5] 
Fig. [2] shows a contour plot of log 10 N cri t- 



3.2 Interpretation of the results 

In order to attempt to interpret the results in the previous 
subsection, it is helpful to illustrate the mass density dis- 
tribution of each component within the cluster and this is 
done in Fig. [3] 
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Figure 3. Illustration of mass distribution in star clusters. The 
dashed line represents r^, the lines at 135 degrees in the cen- 
tre represent the area dominated by the heavy component (i.e 
£4 3> 1), the lines at 45 degrees in the far halo represent the 
area dominated by the light component (i.e B < 1), the crossed 
section represents the area where there is a mixture of heavy and 
light components (i.e. £2- ~ 1 ) 



V/ 



Table 3. Number of heavy stars (N2) at N c rit m units of 10 4 




Figure 4. Left: System with > 1 and with large enough 2^ 
to remove most of the light component from within the half mass 
radius. Right: System with fj 3 - < 1 and ^2 > 1. 



Firstly, let us consider models in which ^ > 1- In a 
region where both components are present at comparable 
densities, there is a strong tendency towards mass segrega- 
tion. Therefore, in the region at which — ~ 1, the ratio — 

' t> pi pi 

is a rapidly decreasing function of the radius, i.e .the transi- 
tion region is narrow. Inside this region, 7712 dominates, and 
mi dominates outside. Clearly the radius at which this re- 

and must be near ru when 
the 



gion is located increases with j^- 

f& = 1 (Fig. [Z]). Finally, for models in which ^ > 1, 
tendency towards mass segregation decreases, the decrease 
of £J with r is more gradual, and the transition region is 
more extensive (Fig. j5J|. For the same reason the regions 
dominated by a single component are more restricted than 
when 2a > 1. 



3.2.1 Dependence on the number of heavy stars N2 

The values of N2 at N cr u are given in Table [3] The varia- 
tion in N2 is considerably less than that of N cr u . For large 
^ and fixed the value of N2 has approximately the 
same value, independent of Now, we give a possible in- 
terpretation of this empirical finding that the stability of the 
system is dominated by the heavy component. 



«2 

Mi 



1.0 


0.57 


0.50 


0.48 


0.47 


0.45 


0.41 


0.35 


0.5 


0.44 


0.40 


0.39 


0.36 


0.34 


0.32 


0.30 


0.4 


0.38 


0.38 


0.35 


0.34 


0.36 


0.29 


0.26 


0.3 


0.34 


0.33 


0.32 


0.31 


0.29 


0.27 


0.25 


0.2 


0.27 


0.28 


0.26 


0.27 


0.24 


0.22 


0.22 


0.1 


0.18 


0.19 


0.21 


0.20 


0.22 


0.18 


0.20 




2 


3 


4 


5 


10 


20 


50 



Firstly, let us consider the case of > 1 and ^ 3> 1 



m 1 



(Fig. [4] left). Within rn the heavy component dominates and 
most of the light component is removed to the outer halo. 
In this case, the light component acts as a container for the 
heavy component. Here, the stellar mass of the light com- 
ponent is not the most important factor, rather the most 
important factor is the overall mass of the container. If we 
were to replace this by an equal mass of stars with stellar 
mass m.2, the behaviour of the stars inside r^ would be nearly 
the same, and so the value of N2 at the stability boundary 
would be roughly the same as for a one-component model. 
Indeed, since part of the container consists of stars of mass 
mi, this could also explain why the values of N2 are in fact 
somewhat less than the value of N cr u for a one-component 
system (i.e. 7000, |Goodma"n| ( |1987[ )) and, in fact, why the 
critical value of N2 is decreasing with decreasing On 
the other hand, the fact that the critical value of N2 is less 
than 7000 may also partly be due to the fact that the energy 
generation rate, equatio n |Tj), is sm aller than that used by 
Goodman. In his paper ( Goodman||1987 equation 11.26) he 
shows implicitly that the critical value of N is approximately 
proportional to the square root of the numerical coefficient 
in e. At any rate, the arguments we have presented are con- 
sistent with the results in the uppermost rows of Table [3] 

Secondly, consider the case < 1 (Fig. [i] right). If 
the system is Spitzer unstable, the heavy component de- 
couples from the light component and forms its own sub- 
system. This heavy subsystem can itself become gravother- 
mally unstable and exhibit a temperature inversion in the 
same way as a one-component model. In this case, however, 
there is not enough mass in the heavy component to domi- 
nate throughout the region within Th, and so it is not quite 
as easy to relate this to the one-component case. Rather, we 
assume that the heavy component behaves like a detached 
one-component model. However, the basic conclusion is still 
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the same, the stability of the model is determined by the 
heavy component. Since the heavy component is again sit- 
ting in the potential well of the lighter stars, it is easier for 
a nearly isothermal region to be set up in the heavier stars 
than if the entire system consisted of heavy stars, and we 
again expect N C ru to correspond to a lower value of JV2 than 
in the one-component case. 

There is also a noticeable increase in the values of -/V2 
with decreasing ^ in the top rows of Table |3| There is 
currently no clear interpretation of this effect Dut it may 
possibly be related to the effect of mass segregation, as the 
region dominated by the heavy component is larger for larger 
% (see Fig. [5). 



Table 4. 



Comparison of values of e and ^ 



m 2 
mi 


M 2 
Mi 




N 


Kim et al log e 


Gas model loge 


2 


0.02 


3 


x 10 4 


-1.620 


-1.553 
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0.03 


3 


x 10 4 


-1.224 


-1.167 


3 


0.03 




10 5 


-1.597 


-1.544 



mi 


M 2 
Mi 




JV 


Kim et al — 

r h 


Gas model — 

r h 
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0.02 
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7.03 x 10~ 3 


4.86 x 10~ 3 
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0.03 
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x 10 4 


1.31 x 10~ 2 


0.91 x 10~ 2 


3 


0.03 




10 5 


5.38 x 10~ 3 


3.93 x 10~ 3 



3.2.2 Goodman's stability parameter 



Goodman ( 1993 1 suggested that the quantity 



Table 5. Value of loge for largest stable run 



(6) 



^ _ Etot/trh 
E c I trc 

should indicate the stability universally, where log 10 e ~ —2 
is the stability limit below which the cluster would become 
unstable. Here Etot is the total energy, E c is the energy of the 
core, trc is the core relaxation time and t r h is the half mass 
relaxation time. |Kim, Lee fc Goodman ( 1998 1 carried out 
research using a Fokker-Planck model which seemed to sup- 
port the condition, although the models they studied were 
all Spitzer stable. 

We have compared the values of e found by |Kim, Lee] 
|fc Goodman| |l998) to results obtained from the gas code 
(Table |4j). All the models co mpared in Table |4} wh ich are 
the same as those studied by Kim, Lee & Goodman ( 1998 1, 



are stable in the post-collapse expansion as well as being 
Spitzer stable. An important difference between the Fokker- 
Planck model used by Kim, Lee & Goodman ( 1998[) and the 
gas code used in this paper is that |Kim, Lee fc Goodman| 
( 1998 1 included an energy generation term in both compo- 



nents, whereas the gas code only contains an energy genera- 
tion term in the heavier component. Therefore, it would be 
expected, in the case of the gas code, that the core would 
have to collapse further in order to generate the required 
amount of energy (from Henon's principle, see Section [3. 3[ ). 
This could explain the differences in the values of — in Ta- 

However, as ^ increases, the heavier component will 
dominate in the core and the energy generation of the lighter 
component will then become negligible. As can be seen in 
Table [4] there is good agreement between the two results 
for log 10 e even though there are only small values of . 



Also, it is possible that Kim, Lee & Goodman ( 1998 1 used 



a different definition of t rc than the one used in this paper 
(see equation [7|. However, as there is such good agreement 
between the values in Table |4j it is unlikely that |Kim, Lee] 
fc Goodman| ( |1998 1 used a significantly different definition. 
Unlike the other models in the present paper, the models in 
Table [4] are Spitzer stable. These runs have only been car- 
ried out in order to make a comparison of the calculation of 
e and ^f- using the gas code with the results of 



Kim, Lee & 



Goodman (19981. Next we will test the use of epsilon as a 
stability criterion for Spitzer unstable cases. 

We tested the stability criterion based on equation (JgJ) 



Ml 
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-2.15 
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for the subset of models given in Table jHJ For each fixed ^ 
and — , the values of e were found to decrease with increased 

mi ' 

N up until the post-collapse evolution became unstable. The 
values of log 10 e given in Table [5] are the values for the run 
with the highest stable N. As can be seen from Table [5] the 
value of log 10 e is indeed in the region of —2. However, the 
limiting value of stable e varies with ^ and to a much lesser 
extent with j^-. 

Now we shall try to improve upon the definition of e. In 
equation pj, t rc and t r h are defined by 



trc — 



G 2 m,2pc,2 In A 



and 



trh. — 



(Grh) 5 In A' 



(7) 



(8) 



Note that t rc was defined using the properties of the heavy 
component in the core rather than the averages of both com- 
ponents, as the heavy component dominates in the core. 
However t r h (equation Q) depends on JV and m, which can 
vary dramatically with and for fixed N2 whereas, as 
argued in Section |3.2[ the important criterion is the number 
of heavy stars. We suggest that a modified version of the 
Goodman stability parameter could be constructed using 
a relaxation time based on the heavy component in place 
of t r h- For example, if |g- > 1 and we assume that the 
heavier component dominates within ru, then the proper- 
ties of the system within rn would be roughly similar to 
that of a one-component system with the same total mass. 
We can attempt to treat the system as if it consisted en- 
tirely of the heavy component with an effective number of 
stars N.t = — . The half mass relaxation time of this one- 
component system would be 
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Table 6. loge and log £2 for the case 



M 2 



loge 


-1.65 
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-2.15 


-2.75 


log £2 


-1.51 


-1.39 


-1.53 


-1.56 
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trh,ef — 



l + 



A/ 2 
Mi 



iU, 



+ 



)( 



In A TV 



™a / \ In A 7V. 



(Gm 2 )5 lnATV e / 

We can define a modified stability condition by replac- 
ing t r h with t r h, e f in the definition of e which would then 
give the following condition 

Etot /trh,, 



£2 = 



Ec/tr 



The values of log 10 e and log 10 £2 are compared in Table[6]for 
the case ^ = 1. The values of log 10 e 2 are in much better 
agreement with each other than those of log 10 e and suggest 
that a better stability condition is log 10 e 2 ~ —1.5 rather 
than log 10 e ~ —2. For the cases with ^ 1 it is unclear 
how to define an appropriate relaxation time, and so we will 
not consider the modified stability condition for those cases. 

To summarise, the values of e (and especially £2) seem 
to give an indication of stability for the two-component mod- 
els but the values of e were found to change with different 
conditions (e.g ^")- The critical value of e 2 is much less 
variable. The critical value of log 10 e or log 10 e 2 is still to be 
tested for multi-component models, and this would be an 
interesting topic for further research. 



3.3 Weak oscillations 



Henon ( 1975 I suggested that the energy generation rate of 
the core is determined by the requirement that it meets the 
energy demands of the rest of the cluster. This demand is 
normally thought of in terms of the energy flux at the half 
mass radius. We shall refer to this as Henon's principle. This 
principle, together with the notion of gravothermal instabil- 
ity, is the basis of the usual qualitative picture of gravother- 



mal oscillations (Bettwieser & Sugimoto 1984), which we 
now recap. 

In a situation with very large TV, the core has to col- 
lapse to a small size in order to meet the required energy 
generation. The steady state is gravothermally unstable, as 
there would be a large density contrast in a nearly isother- 
mal region. If the core is generating more energy than can 
be conducted away, this would cause the core to expand, 
cool and reduce its rate of energy generation. If there is suf- 
ficient expansion, then the core would be cooler than its sur- 
roundings. This would result in the core starting to absorb 
heat. Since the core has a negative specific heat capacity, 
this would cause the core to expand further and became 
even cooler than before |Bettwieser & Sugimoto 1984J) . Ul- 
timately, however, the core must collapse again to meet the 
energy requirements of the rest of the cluster. Here, we adapt 
this explanation of gravothermal oscillations to the case of 
two-component clusters. 



In one-component gas models, as TV increases the insta- 
bility first appears in the form of periodic oscillation:^] ( Heg- 
gie fc Ramamani |1989 1. In order to study the instability for 



the case of weak or low amplitude oscillations in our two- 
component model, a model was chosen which demonstrated 
periodic oscillation with parameters ^ = 2, = 1 and 
TV = 2.0 x 10 4 (the value of TV crit for ^ = 2, |[* = 1 is 
1.7 x 10 3 from Table [2}. Fig.[6]plots In p at various fixed val- 
ues of log r for this model. The total energy flux L is shown 
in Fig.[7]over the particular expansion phase from 24.54ti, r h 
to 25.18^,-h and the contraction phase from 25.18t itr h to 
26.52ti, r h- Fig. [8] shows the profiles of logp and logcr 2 over 
the expansion phase from 24.54t; jr h to 25.18ti, r h- 

During the expansion of the core, the flux in the inner 
region (between r c and rh) drops and eventually becomes 
negative (Fig. [7] top) in a small range of the radius. At this 
point there is an inwards flux of energy to the core. Since the 
core has a negative heat capacity, it would be expected that 
this would enhance the negative flux and therefore the ex- 
pansion. However, the expansion stops at this point. This is 



similar to behaviour observed by McMillan & Engle (19961. 
Now we explain why this happens. Henon ( 1975 1 argues that 
the flux at must be maintained, and we note that there is 
always a positive flux at the half-mass radius rh- Since the 
flux from the core becomes negative at some radius between 
r c and ru, there must be a positive flux gradient in some 
region between the core and half mass radius. This can be 
seen in Fig. [7] (top) towards the end of the expansion and it 
continues into the early part of the contraction phase (Fig. 
[7] bottom). 

The flux gradient can be related to density via equa- 
tion j5J|. As the heavier component dominates in the inner 
regions (see Fig. [8jl , the main contribution to the flux is from 
the heavier component (i.e. L ~ L2). Outside the core the 
energy generation will be negligible. Finally, the temporal 
change in lnp is greater than that in In a 3 . Taking all of this 
into account and rearranging equation |5| will result in the 
following: 



dL 



(- 



hi 



f>2 



(9) 



7^p20"2 

Since all of the coefficients of the flux gradient are posi- 
tive the sign of flux gradient must be the same as that of the 
Lagrangian derivative of the density. Thus a positive radial 
flux gradient in space implies that the density is increasing 
with time. This can be seen in Fig. [6] where the dashed 
lines mark the moment when the contraction becomes an 
expansion, and the solid lines mark the time when contrac- 
tion resumes. It is clear that the contraction begins at large 
radii (logr > —1.6) while the core is expanding, and that 
this region of contraction propagates inwards at later times. 
This can be related to the position of the positive gradient 
in Fig. [7] via the above equation (as long as the density is low 
enough that energy generation is negligible). Therefore, the 
collapse of the parts of the cluster between the core and 
starts while the core is expanding, and brings the expansion 
to a halt. Note that Fig. [6] is density plotted at fixed radius 
whereas time-derivatives in equation [9] are at fixed mass. 

1 Strictly, only periodic if one scales out the steady expansion 



Gravothermal oscillations 7 
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Time {t i rh ) 
Contour plot of log/3 over logr vs Time 




Expansion phase from 24.54?^^ to 25.18ii :r fr 



Time 




logr 



Collapse phase from 25.18fi )r h to 26.57ii )T -/j 




log?- 



Figure 6. Values of log p at fixed values of log r in the range 
-2.7621 to -0.5038 in equal steps of size 0.1129, for two- 
component models with = 2, |^ = 1 and N = 2. Ox 10 4 . Bot- 
tom: contour plot of log p; the dashed lines represent the point of 
highest density reached locally over the time interval 23. 5t^ r h to 
26. 7i; r /, and solid lines are the regions of lowest density reached 
between the dashed lines. X marks the points of core bounce, 
where the core stops contracting and starts expanding 



Nevertheless in Fig. [6] we can also see that there are inter- 
mediate radii in which the density evolves in the opposite 
way from the core. 

Although we have constructed the details of this de- 
scription in the context of two-component models, nothing 
we have said depends entirely on this, and it seems likely 
that similar ideas will apply to one-component and multi- 
component models. 



4 CORE COLLAPSE TIME 

While it may seem that the study of core collapse times is 
inappropriate in the context of gravothermal oscillations, it 
can be argued that the collapse phase of a gravothermal os- 
cillation is not essentially different from the phenomenon of 
core collapse. Furthermore, another reason for its inclusion 
is that the evolution of isolated two-component models is 
an interesting research topic in its own right, and with the 
aim of constructing a comprehensive approximate theory of 
these models, studying the core collapse time is an appro- 
priate first step. 

The core collapse time i cc for a one-component cluster 



Figure 7. Total energy flux in a two-component model with 
HU. = 2, |^ = 1 and N = 2.0 X 10 4 . The expansion phase 
from 24.54ti r /, to 25.184; r /, is on the top and the contraction 
phase from 25.18t; jr /, to 26.57i/ r /, is on the bottom 



with Plummer model initial conditions has been found to be 
approximately 15.5ti jr ./i(Binney & Tremaine: 2008 Heggie & 



[rlut][2003l l using various methods. |Takahashi| ^ 1995[ ) founc 
a longer t cc of 17.6t itr h with a one-component anisotropic 
Fokker-Planck code. However, the presence of a range of stel- 
lar masses can have a dramatic effect on the collapse time be- 
cause of the process of mass segregation. The effect of mass 
segregation in multi-component models has been studied us- 



ing Fokker-Planck calculations (Murphy et al 19901 Chernoff 



& Weinberg 19901 and Monte Carlo methods ( |Gurkan et 



al 2004). The effect of mass segregation in two-component 



models has already been studied extensively using direct N- 
body methods ( jKhalisi et al||2007| ). 

For the gas model runs discussed in Section 3, Table [7] 
gives the values of the collapse time in units of the initial 
half mass relaxation time. Fig. [9] shows a contour plot of log 
t cc . The fastest collapse times occur with models of low ^ 
and high ^. 

For two-component s ystems, the timesca le of mass seg- 
regation varies as (Fregeau et al 2002 and references 
therein) . As mass segregation ennances tne central density, it 
is expected that the mass segregation timescale is compara- 
ble with the timescale of core collapse. Fig. [10] compares the 
variation of the timescale of core collapse with the expected 
timescale of mass segregation. For the case of |j£ = 1.0 (top 
lines in Fig. \W\ the collapse time indeed appears to vary as 



/ TT?2 ^ 

^ mi ' 



However, for lower values of ^ . 

Ml 
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logpi and log/?2 at 24.54^^ and 25.18ti.rZt 




logr 



logo-f and logoj at 24.54t;.,.|, and 25.18ti jr h 




-4 -3 -2 -1 

log?' 



1 



Figure 8. Profiles of log 10 (p) (top) and log 1 g(cr 2 ) (bottom) for 
each component at maximum (dashed line) and minimum (solid 
line) expansion over times shown in Fig. [7] The heavy (light) 
curves refer to the more (less) massive component 



Table 7. Collapse time t cc in units of the initial relaxation time 
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time decreases more quickly than for ( ! 



Khalisi et al 



1 2007) also found a steeper decrease of the core collapse time 

Mi.t, 



0.1, where M^tot is the 



in their study, for the case 
initial total cluster mass. 

We can attempt to improve on these ideas at least qual- 
itatively by considering in more detail a Spitzer unstable 
model. In that case, we can separate the pre-collapse evolu- 
tion of the cluster into an initial mass segregation-dominated 
stage and a later gravothermal collapse-dominated stage, in 
which the centrally concentrated heavy component behaves 
almost as a one-component system thermally detached from 
the lighter component. We propose that this separation can 
be located via the minimum of the rate of change of cen- 
tral density ratio (i.e mm { (^) })■ The reasoning behind 
this is as follows: as time passes, the increase in the den- 
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Figure 9. Contours of log( icc ) as a function of M2/M1 and 



sity ratio caused by mass segregation starts to slow due to 
a combination of decreasing relative density and increasing 
temperature of the lighter component in the central regions. 
We assume that it is at this point that the gravothermal 
collapse of the heavier component becomes the dominant 
behaviour of the system. The gravothermal collapse in the 
heavy component increases the temperature of the heavy 
component, and because the light component absorbs en- 
ergy from the heavy component, the collapse of the heavy 
component causes a deceleration in the collapse of the light 
component. This in turn enhances the rate of increase in the 
density ratio. 

shows the density ratio pz/pi vs time for TV = 
1, 0.1. For the case of ^ = 2 (the lowest 



Fig. 
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10000 and 



~M 2 

Mi ~' ~~ mi 

curve) there is a clear distinction between the part before 
the point of inflection at about j^— = 5 (i.e the initial mass 
segregation phase) and the part after the point of inflection 
(i.e the gravothermal collapse phase). As ^ increases the 
initial phase dominated by mass segregation becomes more 
substantial and eventually the initial mass segregation phase 
brings the system all the way to core bounce. However, as TV 
increases, binary energy generation becomes less efficient rel- 
ative to the energy demands of the cluster ( Goodman|1987 1 . 
Therefore, the core needs to reach a higher density at core 
bounce for larger TV. As the initial phase of mass segrega- 
tion is self limiting for the reason given above, mass segre- 
gation cannot increase the central density beyond a certain 
point. Therefore, it would be expected that the gravother- 
mal collapse dominated phase must eventually return with 
increasing TV for any given j^- and 



5 DIRECT TV-BODY 



Bettwieser & Sugimoto ( 1985 1 compared TV-body systems 
to gaseous models using a direct TV = 1000 model. Even 
though the value of TV is small by today's standards there 
was still fair agreement during the pre-collapse phase. There 
were large statistical fluctuations in the post-collapse phase 
and this was most likely due to the small particle number. 
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Table 8. Collapse time t c 
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Figure 10. Solid lines are log( fcc ) vs log from top to bot- 
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(top) and ^ =0.1 (bottom) Curves from bottom to 
^ = 2,3,4,5, 10,20 and 50 
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Figure 12. log 10 r c vs time, where r c is defined as in NBODY6, 
™a = 2 , Ma = 1 and N = 32fc. 

m i ' M i 



However, it is still important to confirm a sample of the 
results of the gas model by using a direct iV-body code. 
The case of — = 2 4^ = 1 was chosen because it 

mi Mi 

had the smallest value of N cr i t . The values of N used for 
these runs were 8fc, 16fc, 32fc and 64A;. The collapse times of 
the runs in iV-body units (sec Hcggic & Mathieu 1986 ) and 
units of ti : rh are given in Table |8j The average collapse time 
measured in units of U t rh is about 7.5 which is lower than the 
predicted value of 8.95 in Table [7| The difference in collapse 
time could be because of the approximate treatment of two- 
body relaxation in the gas model, the neglect of escape, or 



parameter choices in the gas code (Section 2.1 1. 

For the case of the runs with N equal to 8k and 16k no 
behaviour was found which could be described as gravother- 
mal oscillation. This is in agreement with the gas code, which 
gave N cr it = 17000. However, the 32k case does show a cy- 
cle of expansion and contraction of the core over the time 
interval 4500 to 5500 iV-body units (see Fig.|12[). In order to 
check that the expansion was not driven by sustained binary 
energy generation, we consider the evolution of the relative 
binding energy ^S-, where Ef, is the total binding energy of 
the binaries and E is the absolute value of the total energy 
of the cluster, over this time period. This is plotted in Fig. 
[l3]along with the core radius. There are small changes in the 
binding energy of binaries over this period, decreases as well 
as increases, but this cannot fully account for the expansion 
phase that is observed, as there are other periods with sim- 
ilar binary activity in which no sustained expansion occurs. 
Also, the time scale of the expansion is much longer than the 
relaxation time in the core (~ 0.5 in iV-body time units). 
Therefore, we assume that the expansion must be driven by 
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logr c and relative binding energy vs time 
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Figure 14. 32fc TV-body results, log p in Lagrangian shells of 1, 2, 
5, 10, 20, 30, 40, 50, 62.5, 75 and 90 percent mass in the heavier 
component. It can be seen that the collapse at 5400 starts further 
out (while the core is still expanding) and propagates towards the 
core. 



visible towards the end of this expansion and this is what is 
driving the expansion. From the results of the 32fc and 64k 
runs it seems that the value of N cr u = 17000 obtained by 
the gas code is a reasonable indicator of stability for the TV- 
body case in the sense that none of the signs of gravothermal 
behaviour were found for TV < 16k. 



log^ c 



Figure 13. 32fc TV-body results. Top: relative binding energy 
^ compared to log r c over time 4500 to 5500 TV-body units. 
Bottom: logr c vs logti^ over the same time period, where, • and 
■ represent the starting and finishing points. 



phenomena outside the core, and gravothermal behaviour is 
a plausible explanation. 

Several other pieces of evidence point to this conclusion. 
Fig. |14| shows the density in Lagrangian shells of the heavier 
component. As discussed in Section 3.3 (e.g. Fig.Jfi] top) the 
region further away from the core is seen to contract while 
the core expands. Also, in the cycle of lnp c vs the core ve- 
locity dispersion logt>;?, the temperature is lower during the 
expansion where heat is absorbed and higher during the col- 
lapse where heat is released (Fig. 13 bottom). This is similar 



to the cycles found by Makino 



11996) for one-component 



models and is another sign of gravothermal behaviour. The 
results from the 32k gas run are shown in Fig. [15] for com- 
parison. 

The 64k run shown in Figs. |16|and|17|has large ampli- 
tude oscillations. There is a part of the expansion which is 
shown in Fig. [17] between 7353 and 7390 in which the rela- 
tive binding energy of binaries is nearly constant. Therefore 
bin ary activity cannot be what is driving the expansion. Fig. 
(bottom) shows the evolution of the profile of log v 2 over 
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6 CONCLUSIONS AND DISCUSSION 

The main focus of this paper has been on the gravothermal 
oscillations of two-component systems. The critical value of 
TV for the onset of instability has been found for a range of 
stellar mass ratios and total mass ratios using a gas model. 
The case of = 1 and — = 2 was further investigated 

Mi mi ° 

using the direct TV-body code NBODY6. The value of N cr u 
obtained from the gas code seems to be a good indicator for 
stability in TV-body runs for this case. Based on this, it is 
a reasonable assumption that the other N cr it values would 
give an indication of the stability for direct TV-body systems. 
The values of N crit for the two-component model were found 
to be much higher than for the one-component case and 
were found to vary with ^ and jj 2 -. However, the value of 
TV2 at the stability limit was found to vary much less than 
TV itself. This seems to suggest that instability depends on 



the properties of the heavy component (see 3.2|. A possible 
explanation of this is given in Section [3. 2| 



part of the expansion. A negative temperature gradient is 



The physical manifestation of the oscillations was inves- 
tigated for the case of small-amplitude periodic oscillations 
in the gas model. It has been pointed out that the collapse 
of the region between r c and is an important mecha- 
nism which can halt the expansion phase of a gravothermal 
oscillation. This mechanism should also be present in one- 
component models and it would be an interesting topic for 
future work to see how this mechanism would behave with 
different stellar mass functions. 

|Kim, Lee &: Goodman] |l998) argued that two- 
component clusters may be realistic approximations of 
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Figure 15. Evolution of the central properties of the heavy com- 
ponent in a 32k gas model. Top: logp2,c vs time (units t; jT -fi), 
bottom: logp2, c vs logu| c . All cycles are clockwise. The initial 
drop in log v\ results from the two components trying to achieve 
thermal equilibrium 

multi-component clusters, where the two components are 
neutron stars and main sequence stars and the effect of 
white dwarfs (heavier than the turnoff mass) was assumed 
to be negligible. They also only studied cases that were 
Spitzer stable, which means that the components were able 
to achieve equipartition of kinetic energy. For the two- 
component case, it is only possible for it to be Spitzer stable 
if there is only a small amount of the heavier component 
present. As there is a significant range of stellar masses in 
a real star cluster, it is likely that some form of the Spitzer 
instability will be present. 

To apply our ideas to a multi-component system, it may 
be possible to group the heavier components together if they 
are able to achieve approximate thermal equilibrium. This 
could be considered as a single heavier component which is 
Spitzer unstable with respect to the remaining components. 
This would help to reduce a multi-component system to the 
two-component case studied in this paper. 

Nevertheless, it is not clear quantitatively how the con- 
siderations of this research are to be applied to a multi- 
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Figure 16. 64k iV-body run. Top: logp2 at various Lagrangian 
shells in the heavier component, bottom: logr c vs time 

component cluster. Furthermore, we have ignored many 
things such as primordial binaries, tidal fields and stellar 
evolution and these are important in the evolution of a real 
star cluster. Further study is needed in order to understand 
the phenomenon that is gravothermal oscillation. 
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Figure 17. Top: evolution of the relative binding energy ^ and 
logr c . Bottom: logi>| measured at Lagrangian shells in the heav- 
ier component over part of the expansion of the core. 



Figure 18. 64k gas model. Top: logp2,c vs time (units ti. r h), 
bottom: logn| vs logr at different times during a gravothermal 
expansion. 
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